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1. Introduction 


Cardiac arrhythmias represent a leading cause of morbidity and mortality in developed 
countries. In spite of intense research, the mechanisms of generation, maintenance and 
termination of theses arrhythmias are still not clearly understood. The ability to predict, 
prevent and treat these arrhythmias remains a major scientific challenge. In the last years, 
new technologies have produced a huge amount of information at different levels: 
subcellular, cellular, tissue, organ and system. However, the knowledge of each 
component is not sufficient by itself to understand its behavior and an integrative 
approach, which should take into account the complex interactions between all 
components is needed. This understanding is necessary to improve prevention and 
treatment of cardiac arrhythmias. 

Traditionally, cardiac research has been based on experimental and clinical studies. 
However, in the last decades, a quantitative understanding of the biophysical and 
biochemical properties of individual cardiac myocytes and of the anatomical structure of the 
heart has permitted the development of computational models of single myocytes, portions 
of tissue and even the whole heart. These models have opened a new approach in studying 
the complex mechanisms underlying cardiac arrhythmias by means of computer 
simulations. 

Integrated and multiscale models of the heart are improving very fast, as a large body of 
research is being made in this field. These models integrate, with a great degree of detail, the 
different levels of the organ, from the genetic characteristics of membrane ionic channels, to 
the electrophysiological characteristics of cardiac cells and to the anatomical structure of the 
different cardiac tissues. In this chapter, the basis of cardiac multiscale modeling will be 
reviewed and recent and specific case studies for the different levels will be described 
focusing on the understanding of cardiac arrhythmias during myocardial ischemia and drug 
treatment. 


2. Modeling cardiac ion channels activity 


Comprehensive models of the electrical activity of cardiac ion channels were formulated for 
the first time in 1952 (Hodgkin & Huxley, 1952). These models were based on the Hodgkin- 
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Huxley mathematical formalism, which allows the reproduction of macroscopic ionic 
currents. This approach computes the conductance of a certain ion channel as a function of 
the open probability of each gate of the channel and the maximum conductance of the 
membrane for that particular ion. In additon, gate transitions from closed to open and vice 
versa follow a first order voltage dependent behavior at a rate that is independent of the 
remaining gate states. 

The improvements in the experimental techniques employed to measure and characterize 
the ionic currents responsible for the cardiac action potential have provoked a paralell 
evolution in the ion channel models, which have acquired an increasing complexity and 
realism. Two examples are the separation of certain currents into its components, like the 
slow and rapid delayed rectifier potassimum currents, and the discovery of ion current 
dependencies on electrolyte concentrations. Importantly, experimental data have evidenced 
the existence of state dependent transitions of ion channels that Hodgkin-Huxley formalism 
fails to reproduce. This fact has favored the development of Markov models for ion 
channels, which are based on the assumption that the transitions between channel states 
depend on the present conformation of the channel and not on its past history. Furthermore, 
recent developments in molecular biology and in the genetics of ion channels have made 
possible the characterization of the electrophysiological properties of a certain mutated ion 
channel (Tomaselli et al., 1995), allowing the development of its mathematical model using 
Markov models. Although Markov models can provide a more realistic formulation of ion 
channels, Hodgkin-Huxley models are also currently used because their computational 
requirements are significantly less demanding. Both formalisms are used to model ion 
channel-drug interactions and pathological situations. 


2.1 Modeling myocardial ischemia at the ion channel level 

The ionic current and action potential models developed in the last decades make it possible 
to simulate the electrical activity of cardiac cells both in normal and abnormal conditions. 
Among the latter, acute myocardial ischemia is one of the most threatening abnormalities 
that can occur in cardiac tissue because it can easily trigger malignant and potentially mortal 
reentrant tachyarrhythmias such as ventricular tachycardia (VT) and ventricular fibrillation 
(VF) (Janse & Wit, 1989). Indeed, VF subsequent to ischemia is considered as the first cause 
of mortality in Europe, the USA and a significant part of Asia (Ross, 1999). 

Myocardial ischemia occurs when a coronary artery is occluded. When this happens, a 
certain zone of the myocardium (formed by the cells that depended on the occluded artery 
for their blood and nutrient supply) begins to suffer metabolic and electrophysiological 
changes. Specifically, K+ ions begin to rapidly accumulate in the extracellular medium (a 
phenomemon termed hyperkalemia), the media (both intra and extracellular) become 
acidotic (with a drop in pH from 7.2-7.4 to 6.2-6.4 in 10-15 minutes), and oxygen levels 
decline (hypoxia). Hypoxia provokes in turn a depletion in intracellular ATP levels and a 
concomitant increase in ADP levels, which results in the partial activation of a specific type 
of K* current called ATP-sensitive potassium current (Ikatp)) (Noma, 1983). 

To properly and predictively simulate the cardiac arrhythmias arising from myocardial 
ischemia, the model must be of a multi-scale nature and ischemia has to be modeled in the 
ion channel level, the cell level and the tissue level. While tissue modeling of ischemia will 
be adressed further in this chapter, ion channel ischemia modeling should be done as 
follows. First, hypoxia must be included in the model through its effects in Ikiarp) current 
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activation (see below). Second, hyperkalemia must be simulated by increasing [K+]o levels 
(Rodriguez et al., 2002). Elevated [K+], in turn alters the conductance of a number of 
sarcolemmal channels, such as Iki or Ik, and these effects are included in most ion channel 
models published to date. Third, acidosis must be modeled through its effects in the 
maximum conductances and kinetic properties of the Na+ and Ca% inward currents 
(Ferrero, Jr. et al., 2003; Shaw & Rudy, 1997a). 

The Ikiatp) current, which is normally dormant in normoxic situations but is partially 
(though only slightly) activated in acute myocardial ischemia, is carried by K* ions (Noma, 
1983), shows voltage-dependent inward rectification (but is otherwise non-voltage 
dependent) caused by Mg?+ and Nat ions (Horie et al., 1987), and is gated by intracellular 
ATP with intracellular ADP playing a modulating role in the ATP-mediated deactivation 
(Weiss et al., 1992). All these facts have to be present in an ion channel model of the current. 
In 1995, Ferrero et al. published the first comprehensive model for the myocardial Ik(arp) 
current (Ferrero, Jr. et al., 1996), which has been further improved more recently (Michailova 
et al., 2005). 

The mathematical model for the Ikiatp) current includes all the aforementioned effects, 
which are represented in the following equation: 


Ikar) = O 8o Po Vir f arp (Vm- Ex) (1) 


where o is the channel density, go the unitary conductance (which is [K*], dependent), yir is a 
term that includes the inward rectification properties of the channels, farp is a term that 
expresses the ATP and ADP dependence of the current, Vm is membrane potential and Ex 
the K+ equilibrium (Nernst) potential. The inward rectification term is expressed as follows: 


= 1 1 
ir 24 2 [Na+]; (2) 
Ki, Mg h,Na 


where Knog and Knog are voltage and [K+]o dependent. Finally, the ATP-ADP dependence 
follows another Hill-type equation as follows: 


1 


fare = : | FATE wj ©) 
K arp (lADP]) 


where both the semi-inhibition constant Karp and the Hill exponent H depend on the 
intracellular ADP concentration. 


2.2 Modeling drug-ion channel interaction 

Also at the ionic channel level, the effect of antiarrhythmic drugs can be mathematically 
formulated. Models that simulate effects of drugs and reproduce the associated clinical 
electrocardiography signal are highly relevant to assess and predict antiarrhythmic actions, 
and help to improve the efficacy and safety of pharmacologic therapy. Additionally, it starts 
to be considered as an important technological tool for helping the pharmaceutical industry 
to develop new antiarrhythmic drugs spending less time and money. In this field, Dr. 
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Starmer is one of the pioneers in modeling drug effects on sodium channels and an 
increasing number of research groups are developing 3D models of the heart to test drug 
actions. In the last years, our group has carried out different studies of antiarrhythmic and 
proarrhytmic effects of different drugs, such as pinacidil, dofetilide and lidocaine on cardiac 
arrhythmias. 

To simulate the effects of drugs on ion channels, different kinds of models can be used 
(Brennan et al., 2009). Sometimes it is sufficient to change the conductance of ion channels 
by introducing the effect of the drug, but to be able to reproduce use- and frequency- 
dependent block it is necessary to use a state-dependent block model. In the state-dependent 
block model, drug interaction with ionic channels depends on the state of the channel. Drug- 
associated channels (blocked channels) do not conduct. The fraction of blocked non- 
conducting channels depends on the equilibrium between blocked and unblocked states. 
Two hypotheses have been proposed to reproduce the state-dependent block: the 
modulated receptor hypothesis and the guarded receptor hypothesis. The Modulated 
Receptor (MR) hypothesis (Hille, 1977; Hondeghem & Katzung, 1977) states that the affinity 
of a drug changes with the state of the channel. Additionally, the rate constants of state 
transitions are different for bound and unbound satates. This is the most general model to 
reproduce the effects of drugs; however a great number of parameters have to be estimated. 
The Guarded Receptor (GR) hypothesis (Starmer et al., 1984) assumes that the drug has a 
fixed and state-independent affinity for their receptor sites, but its access to the binding sites 
is controlled by the voltage-dependent channel gates. In this hypothesis, the rate constants 
for the state transitions do not change when the drug binds to the channel. The GR model is 
less complex than the MR model for the same drug-ion channel interaction as fewer 
parameters have to be estimated. 

In the present section three different models of drugs: pinacidil, lidocaine, and dofetilide are 
described. 


2.2.1 Modeling pinacidil effects 

Pinacidil is a vasodilator drug and is known to increase Ix(atp) in cardiomyocytes, as well as 
in vascular smooth muscle and pancreatic Beta cells (Fan et al., 1990; Nakayama et al., 1990). 
The activation of ATP-sensitive potassium (Karp) channels in heart tissue, especially during 
myocardial ischemia, markedly changes action potential (AP) configuration by reducing AP 
duration (APD) (Nichols et al., 1991) and plays an important role in cardiac 
arrhythmogenesis (Janse & Kleber, 1981). To analyze the controversial effects of this drug 
through computer simulations, the first step is to build a model of the effects of pinacidil at 
the ion channel level. Pinacidil enhances the activity of Ik(atp) increasing the fraction of open 
Karp channels (farp) (Fan et al., 1990). The equation of Ik(arp) formulated by Ferrero et al. 
(Ferrero, Jr. et al., 1996) was modified, as indicated in the following equations: 


Ikat = © 8o Po f arp fein (Vin- Ex ) (4) 
P H 
fare =1 Pick a l Tika [P]+ k [PID] (5) 
[Al? + 1%11°13 11°13 1°13 1 
ky kak + kikis [D] ie k,k3[P] + k,[P][D] 
kay kao k43 kaa 
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The equation of farp was based on the pharmacokinetic scheme depicted in Fig. 1. This 
scheme describes the mechanism of action by which ATP molecules close Karp channels, 
while pinacidil molecules open it. ATP, ADP and pinacidil molecules bind to different sites 
in the protein receptor (Nichols et al., 1991; Weiss et al., 1992). According to the open states 
when molecules of pinacidil and/or ADP are bound to the protein channel, the fraction of 
activated Karp channels was expressed in equation (5). Experimental data (Fan et al., 1990; 
Nakayama et al., 1990; Weiss et al., 1992) were fitted to equation (5) to determine the values 
of equilibrium dissociation constants and the Hill coefficient (H) (see Trenor et al. 2005 for 
details). A factor (fpr) addressing the voltage-dependent block of Karp channels by pinacidil 
was also included in Ix(arp) equation, as in Fan et al. (Fan et al., 1990). 


_ +P o 
kir 
+HA +P +HA 
RAu = R — PR - PRA 
ka ky kas 
+D k +D kı +D k +D ky’ 
+HA +P HE 
DRAs <~ =- DR + PDR => PDRAs 
l d2 12 d4 
+P 
kız’ 


Fig. 1. Cyclic scheme representing different possibilities of binding to the channel receptor 
protein. R, P, A and D stand for the unoccupied receptor, pinacidil, ATP and ADP 
respectively. ki are the overall equilibrium dissociation constants for the different steps of 
the cyclic reaction, and H is the Hill coefficient. 


2.2.2 Modeling lidocaine effects 

Lidocaine is a class I antiarrhythmic drug that exerts its effect by blocking inward sodium 
current (Ina) in a use-dependent manner, and is more effective for high stimulation 
frequencies (Clarkson et al., 1988). The block of lidocaine exists both as neutral and charged 
forms at a physiological pH. The charged form predominates at low pH (6.4) due to the 
protonation of the neutral form with hydrogen molecules. Experimental data have 
confirmed that lidocaine is more effective and recovery from block is slowed (Broughton et 
al., 1984; Moorman et al., 1986; Wendt et al., 1993) under acidosis. 

The formulated model of lidocaine-INna interaction takes into account the modulatory effect 
of pH on the Nat channel block and is based on the GR theory proposed by Starmer 
(Starmer & Courtney, 1986). Three main processes have been considered: the hydrophobic 
pathway, the hydrophilic pathway, and coupling between blocked channels by charged and 
neutral forms using a proton exchange process. Several experimental studies have 
suggested that the block of neutral drug takes place during the activated (O), inactivated (I) 
and resting (R) states (Bean et al., 1983; Clarkson et al., 1988), whereas the interaction 
between the charged form with the channel only takes place in the activated state (Liu et al., 
2003; Moorman et al., 1986; Yeh & Tanguy, 1985). Fig. 2 represents the kinetic diagram of the 
complete model, and equations (6) and (7) mathematically describe the interaction of the 
neutral and charged form with Na* channels, respectively. 
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dt | mh ko +m? (1-hj)-kı +(1-m°)-kr |-[Dy]-(1-by -be) = 


(1) 
| mhj Io +m? (1=hj)- 1, +(1-m)-Ig |-by +l be —ky[ H* |-by 


db aE 
eo Lehi ke |:[Dc]-(1- by ~ be) = neni ft ate be—Ip tet hl] (2) 


by and bc stands for the fraction of channels blocked by the neutral and charged form. The 
association and dissociation rate constants for the different states (activated, inactivated, and 
resting) are ko, kr, kr, lo, h, and lr, respectively. The activation channel gate is m and the 
inactivaction channel gates are h and j. Additionally, we have taken into account the rate 
constants kp and lp to determine the proton exchange process. In equation (7), the 
association and dissociation rate constants are kc and lc, respectively. 


Neutral | Charged 


z > [Dyl* R 


hi R I] 


i 
NJ [Dy] + O + [Da 


al, ea ei 


by 1-b,-b, b by 


Fig. 2. Lidocaine interaction with Na* channels in both forms of the drug: neutral (bn) and 
charged (bc). 


As stated in the GR hypothesis, we introduced the factor (1-b), where b is the sum of the by 
and bcin the Ina equation. 


2.2.3 Modeling dofetilide effects 

Dofetilide is a specific and potent blocker of the rapid component of the delayed rectifier K* 
current (Ikr) with an ICso in the nanomolar range (3.9-31 nmol/L) for ventricular myocytes 
(Carmeliet, 1992; Jurkiewicz & Sanguinetti, 1993; Weerapura et al., 2002b). Dofetilide is 
classified as a pure class III antiarrhythmic agent and provokes a prolongation of action 
potential duration (APD) without any effect on the resting membrane potential, action 
potential amplitude or maximum rate of depolarization (Jurkiewicz & Sanguinetti, 1993; 
Tande et al., 1990). 

To model the effects of dofetilide on Ikr, the formulation of Ik, proposed by Zeng and Rudy 
(Zeng et al., 1995) for guinea pig ventricular cells was used. This formulation is based on 
experimental data obtained by Sanguinetti and Jurkiewicz (Sanguinetti & Jurkiewicz, 1990). 
In their model, Ikr channels present three possible states: Closed (C), Open (O) and 
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Inactivated (I). The formulation of Ikr: in the ventricular cell model (Zeng et al., 1995) is 
expressed in equation (8). 


Ixy = Cig nae Vx T Ex,) (3) 


where Vm is the membrane potential, Ex, is the reversal potential, Gxrmax is the maximum 
conductance of Ikr Xr is the activation gate and R is the time-independent inactivation gate. 
In our dofetilide-Ix, model, the effect of dofetilide is represented by introducing the factor 
(1-b) in the Ik: formulation (where b is the fraction of channels blocked by the drug). Thus, 
the new formulation of the rapid component of the delayed rectifier K+ current taking into 
account the effect of dofetilide is 


Ikr = (1 ~ b) GxrmaxX+ (Vn z Eķ,) (4) 


The blocking activity (b) on Ik, equation is modeled using the GR hypothesis (Starmer et al., 
1984). Dofetilide experimental studies have shown drug-receptor interaction in open and 
inactivated states but not in closed states (Weerapura et al., 2002a; Yang et al., 1997). It has 
been also observed that the drug is trapped when the channel closes, thereby the transition 
between open and closed state is possible but dissociation of the drug does not occur at 
closed states (Carmeliet, 1992; Weerapura et al., 2002a). Taking into account these 
experimental evidences we suggest the model of Ikr block by dofetilide shown in Fig. 3. 
Using this model, the blocking factor b for a concentration of dofetilide D can be calculated 
as: 


aD OIX Rot X O- RKIDO-b)[X Roi X,(1-R)r Jb (5) 


dt 


— 


Fig. 3. Schematic model of dofetilide interaction with Ik, channels in the different states: 
Closed (C), Open (O) and Inactivated (I; and Iz) potassium channels. k and r are the 
association and dissociation constants for the different states. 


3. Modeling cardiac electrical cellular activity 


The electrical activity of the cell is determined by the interaction of the multiple ionic 
channels and ion exchangers present in the cellular membrane. Their activity conditions 
the cardiac AP, which is intimately synchronized with the mechanical activity of the heart. 
Several computational models of cardiac cellular electrical activity have been developed, 
and characterize the activity of ventricular myocytes, atrial myocytes, sinus node and 
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Purkinje fibers for different animal species. The degree of complexity of these models 
increases going from heuristic models describing the physical process of cell excitation 
and propagation, such as Fitz-Hugh-Nagumo model (Fitzhugh, 1961), to 
phenomenological models (Bueno (Bueno-Orovio et al., 2008) model), and to extremely 
detailed ionic models. The most detailed ionic models for AP include transmembrane 
currents, ion transfer through transporters and dynamic changes of ion concentrations 
during the action potential. The ventricular AP models most extensively used in 
simulations because of their detailed and realistic formulation are Luo and Rudy model 
(Faber & Rudy, 2000; Luo & Rudy, 1994; Zeng et al., 1995) for guinea pig, Shannon 
(Shannon et al., 2005) and Mahajan (Mahajan et al., 2008) models for rabbit, and ten 
Tusscher (ten Tusscher & Panfilov, 2006) and Grandi (Grandi et al., 2010) models for 
human. As regards atrial myocytes, Nygren (Nygren et al, 1998) model and 
Courtemanche (Courtemanche et al., 1998) model have also had an important 
repercussion. Luo Rudy AP model was one of the first detailed models formulated and 
has been a fundamental base on which other models have been built. This model has been 
improved and extensively used to simulate the ventricular electrical activity. In these 
detailed AP models, the formulation of the ionic currents follow Hodgkin-Huxley, 
Goldmann-Hodgkin-Katz (GHK) or Markov formalisms, and the membrane potencial Vm 
is related to the ionic currents through equation (11). 


dV, 
I, = nge Ast Dy (6) 


where the total ionic current through the cell membrane (Im) is expressed as the sum of all 
ionic currents (I, through ionic channels and I, through pumps and exchangers) plus the 
current through the membrane equivalent capacitor. 


Fig. 4 shows a scheme of the main ionic currents of a ventricular myocyte (panel A) and the 
time course of membrane potential and some ionic currents (panel B) obtained by solving the 
non-linear differential equations system defined by equation (11) and the different ionic 
current equations in Luo and Rudy model (Faber & Rudy, 2000; Zeng et al., 1995). 


V 3 pA/cm? 


1 uA/cm? 
l 
CaL 


Cleft Space 


-1 Wace 


100 ms 


Extracellular Space -350 A/cm?! Ina 


Fig. 4. Panel A: Schematic representation of a ventricular myocyte and the most significant 
ionic currents through the sarcolemma. Panel B: Temporal evolution of the membrane 
potential Vm, Ca2+, Na* and K+ currents, and currents through the Nat+/K+ pump (Inax) and the 
Nat/Ca?+ exchanger. 


www.intechopen.com 


Multiscale Modeling of Myocardial Electrical Activity: From Cell to Organ 345 


3.1 Modeling myocardial ischemia at cellular level 

Acutely ischemic action potentials can be simulated using these AP models appropriately 
modified to represent ischemia-activated currents (such as Ix(arp)) and ischemic-modified 
currents (e.g. the inward Nat and Ca% currents), as explained previously. 

One of the main components of acute ischemia is hypoxia (i.e. lack of oxygen in the tissue), 
and the main change exerted by hypoxia to myocardial action potentials is a significant 
shortening in APD (Morena et al., 1980). The model proposed by Ferrero et al. for the ATP- 
sensitive K* current (Ferrero, Jr. et al., 1996) yielded the first theoretical proof of the 
determinant role of the Ixiarr) current in hypoxic action potential shortening (Ferrero, Jr. et 
al., 1996), which is an arrhythmic factor due to the APD dispersion it generates across the 
myocardium (Janse & Wit, 1989). Indeed, as shown in Fig. 5, although very few Karp 
channels do activate during acute myocardial ischemia (<1% in the first 15minutes), this 
degree of activation is enough to account for the >50% shortening in APD which is observed 
in experimental conditions (Morena et al., 1980). 

When the simulated myocyte is subject to hyperkalemia and acidosis, as well as hypoxia, the 
main features of the ischemic AP waveforms are correctly reproduced (Ferrero, Jr. et al., 
2003; Shaw & Rudy, 1997a). A shown in Fig. 6, ischemia produces APD shortening, resting 
potential elevation, upstroke phase division into two phases (the first one mediated by 
inward Na* current, the second one by inward Ca% current through the L-type channels). 
Moreover, acute ischemia changes important intrinsic AP properties which are pivotal in 
arrhythmogenesis (e.g., postrepolarization refractoriness). Overall, this ischemic AP 
simulations are the first level of the multi-scale ischemia simulations that will be described 
further in this chapter. 


-100 - 100 ms 


Fig. 5. Hypoxic action potentials simulated using a modified version of the dynamic Luo 
& Rudy model (Luo & Rudy, 1994) in which the Ikarp) model formulated by Ferrero was 
included (Ferrero, Jr. et al., 1996). Numbers indicate percentage of Karp channel 
activation. 


3.2 Effects of drugs at cellular level 

As mentioned before, the electrical activity of the cell is altered under the effects of drugs. 
Once the effects of the drug have been formulated at the ionic channel level, the 
alterations of the electrical activity of the cell can be analyzed including the modified 
current in the AP model. In this way, the effects of pinacidil and lidocaine on ischemic 
cells are described. Also the effects of dofetilide on different types of ventricular cells will 
be analyzed. 
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Fig. 6. Action potentials and inward ionic currents during control conditions (panel A) and 
10 minutes of acute ischemia (panel B). The lower panels correspond to a zoomed view of 
the first 30 ms of the action potentials shown on the top panels. A significant decrease in 
inward Nat current, together with an increase in Ca*? current, is observed. The low Ina in 
ischemia provoques the division of the upstroke phase of the action potential into two 
distinct phases, as seen experimentaly (Janse & Wit, 1989). 


3.2.1 Effects of pinacidil on AP 

As mentioned in section 2, pinacidil shortens the APD by the activation of Karp channels. 
To validate pinacidil model, the formulation of the effect of pinacidil (equations (4) and 
(5)) on Ikare) was introduced into the Luo and Rudy model (Luo & Rudy, 1994; Zeng et 
al., 1995). Simulations were carried out combining different concentrations of pinacidil 
with various [ATP]; and [ADP]; values corresponding to several stages of early acute 
ischemia (Cole et al., 1991; Nakayama et al., 1990; Weiss et al., 1992). The resulting APs are 
shown in Fig. 7. Consistent with experimental findings (Cole et al., 1991; Nakayama et al., 
1990), our results show that increasing concentrations of pinacidil and/or decreasing 
intracellular ATP levels reduce APD. It is noteworthy that changes concerning 
intracellular ATP and ADP concentrations, occurring during acute ischemia, significantly 
enhance the effects of pinacidil. These changes in APD (and hence in cell refractoriness) 
could be key in determining the reentry probability in acute ischemia after the 
administration of the drug. 
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Fig. 7. Set of simulated single cell action potentials in the presence of different 
concentrations of pinacidil: 0, 5 and 10 umol/L from left to right. Shortening of APD can be 
also observed under different concentrations of ATP: 6.8 mmol/L (solid line), 5 mmol/L 
(dashed line) and 2 mmol/L (discontinuous line). 


3.2.2 Effects of lidocaine on AP 

The block of Ina by lidocaine leads to a decrease in the maximal AP upstroke (dV/dtmax). To 
simulate the effects of lidocaine on AP features, we incorporated the model of lidocaine into 
the Luo and Rudy AP model (Luo & Rudy, 1994). Specifically, we measured the reduction of 
the maximum AP upstroke for different stimulation frequencies and pH values when 
increasing concentrations of the drug. We examined the response of dV/dtmax to changes in 
pH from 7.4 to 6.4 for different concentrations of the drug (50 and 100 of umol/L lidocaine) 
and at different stimulation frequencies. Panel A in Fig. 8, shows a decrease of 66% and 83% 
in dV/dtmax from the control value for 50 and 100 umol/L lidocaine, respectively, at a basic 
cycle length (BCL) of 500 ms and pH 6.4. For a pH of 7.4 (see Fig. 8 panel B) , dV/dtmax was 
reduced in 45% and 61% for the same concentrations and BCL values. These data support 
the fact that lidocaine has a greater effect under lower pH. Furthermore, the use-dependent 
effect is more pronounced at lower pH, as the differences in the decrease of dV/dtmax for 
different stimulation frequencies are higher for low pH. 


H 6.4 i 

Peis : p Bi pH 7.4 

r ; | | “Sex. 1000 ms 

Zz a a | je Ra 

ae 0 6 t. oa =n | TNE Zr TS 

g ` Wo o woo oo Ve 

5 | 

ee | cot 

S025 jj Stig presses. 


0 40 80 0 40 80 
Lidocaine (uM) Lidocaine (uM) 


Fig. 8. The effects of lidocaine on dV/dtmax at pH values of 6.4 and 7.4 at different BCLs. 
The steady-state dV/dtmax was normalized (Nor) with respect to its value in the absence of 
the drug. 
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3.2.3 Effects of dofetilide on AP 

The steady-state effects of dofetilide concentration on AP characteristics for the different kinds 
of guinea pig ventricular cells at different BCL values are shown in Fig. 9. In all simulations, 
once dofetilide was applied the cellular models were stimulated until steady- state was 
reached. Only magnitudes associated to the last stimulation pulse are shown. Fig. 9A shows 
that the prolongation of APDoo, for all tested dofetilide concentrations, significantly depends 
on the type of cell. APDo of M-cells was markedly prolonged compared to endocardial and 
epicardial cells. For a concentration of 10 nmol/L, APDoo of endocardial cells was prolonged in 
8.8%, 11.5% and 12.4% of control value (APDoo before the application of drug) for BCL of 300, 
1000 and 2000 ms, respectively. For the same conditions, APDoo of epicardial cells was 
increased in 5.8%, 6.9% and 7.1% of control, respectively. However, under these conditions, 
APDo of M-cells was substantially increased in 14.8%, 36.1% and 292.7% of control, 
respectively. The three types of cells present reverse use-dependency. Due to this effect, the 
prolongation of APDoo increased whith the BCL. Fig. 9B shows the relative increase of APDoo 
(AAPDop in%) induced by different concentrations of dofetilide at different BCLs. 

Fig. 9B shows that dofetilide presents reverse use-dependence in the range of tested 
concentrations (10 nmol/L to 1 pmol/L) with a much more pronounced effect in M-cells 
than endocardial or epicardial cells. Even more, when the concentration of dofetilide 
increased, the reverse use-dependency was enhanced. 
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Fig. 9. Effects of different concentrations of dofetilide on steady-state AP waveform (panel A) 
and APD. prolongation (panel B) from guinea pig endocardial, epicardial and M-cells at 
different BCL values. 
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4. Modeling cardiac tissue electrical activity 


Abnormalities in impulse propagation are involved in the genesis and maintenance of 
cardiac arrhythmias and many experimental and theorical studies have analyzed action 
potential propagation both in normal and pathological situations. The electrical coupling 
between cells and the propagation process is described in the bidomain model (Tung, L., 
1978). This model consists of two continuous domains, the intracellular and the extracellular 
domain, where electrical currents are governed by Ohm's law. The following equations 
describe the bidomain model: 


V-(DVV,,)+V “(DV V,) =Cin i +I), +L (7) 
V -(DVV,,)+ V -((D; +D.)VV,)=0 (8) 


where Dj and De are the intracellular and extracellular conductivity tensors, respectively, Ve 
and Vm stand for the extracellular and the membrane potentials, respectively, Cm is 
membrane capacitance, lion represents the total ionic current density computed using AP 
models, and I, is the stimulus current density. 

Considering that anisotropy similarly affects the intracellular and the extracellular 
conductivity tensors (Di=ADe, à being a scalar variable), the bidomain model can be 
simplified into the monodomain model and the previous equations can be reduced to the 
following one: 


V-(DVVin) = Con Mn Plon dy (9) 
A 
where D= i (10) 
1+4 


Subsequently, the monodomain model consists of an elliptic partial differential equation 
and a parabolic partial differential equation coupled to a system nonlinear ordinary 
differential equations describing the ionic current through the cellular membrane. As this 
model is mathematically simpler and less computational demanding than the bidomain 
model, it is widely used for cardiac electrophysiology simulations. 

Heterogeneous structures of cardiac tissue can also be considered in tissue models allowing 
a wide use of cardiac models in the study of ionic mechanims of arrhythmias. 


4.1 Tissue model of myocardial ischemia 

In the case of pathological tissues, not only the altered AP models have to be considered but 
also the spatial inhomogeneities in the electrophysiological changes distributed within the 
tissue. If we consider the ischemic pathology, after coronary occlusion, the lack of oxygen 
and blood flow provokes important and heterogeneous electrophysiological changes in the 
affected ventricular cells defining an ischemic zone and a border zone, setting the stage for 
reentrant arrhythmias and VF. A large body of clinical and experimental studies try to 
unravel the mechanisms by which acute myocardial ischemia provokes life-threatening 
arrhythmias. However, the rapid changes arising during acute myocardial ischemia and the 
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limitations in experimental techniques hamper the complete understanding of this 
pathology. In this way, computer simulations overcome these limitations providing a 
helpful tool to understand the electrophysiological changes arising during myocardial 
ischemia and the underlying mechanisms of reentrant arrhythmias observed during its 
acute phase. At the tissue level, our group has developed a regional model of ischemia, and 
the consideration of the temporal evolution of ischemia has allowed undertaking 
computational studies to assess the changes in the vulnerability of the ischemic tissue to 
reentry. Not only ionic currents have been analyzed to understand the process of reentry 
initiation but also the traditional indicator for reentry, i.e. dispersion of repolarization, and 
the safety factor (SF) for conduction (a parameter that quantifies the source-sink relationship 
of propagation) as an indicator for conduction block, have been evaluated. 

A regional phase Ia ischemic 2D anisotropic monodomain tissue (Ferrero, Jr. et al., 2003; 
Romero et al., 2009; Trenor et al., 2007) was modeled with a high degree of 
electrophysiological detail. Changes in electrophysiological parameters were modeled 
according to experimental results (Coronel, 1994; Ferrero, Jr. et al., 1996; Irisawa & Sato, 
1986; Weiss et al., 1992; Yatani et al., 1984). These changes are described in Fig. 10A, which 
shows the time-course of the three main components of acute ischemia: hypoxia, 
hyperkalemia and acidosis. Fig. 10B represents the distribution of a simulated cardiac 2D 
tissue comprising a normal zone (NZ), a circular central ischemic zone (CZ), and a ring- 
shaped ischemic border zone (BZ) surrounding the CZ. In the BZ all parameters affected by 
ischemia follow linear spatial gradients, as experimentally recorded (Coronel et al., 1988). 
Membrane kinetics was simulated by using a modified version of the Luo-Rudy dynamic 
AP model (LRd00) (Faber & Rudy, 2000) (see Romero et al. 2009 for details), which includes 
the mathematical formulation of Ixiarr) proposed by Ferrero et al. 1996 (Ferrero, Jr. et al., 
1996). 

A basic stimulus (S1) followed by a prematures stimulus (S2) at different coupling intervals 
(CIs) was applied on the edge of the tissue to quantify the vulnerable window (VW), an 
indicator of the vulnerability to reentry. The VW was defined as the interval of CIs that led 
to reentry. Fig. 10C shows gray-coded voltage snapshots of a figure-of-eight reentry 
developed after premature stimulation at minute 8.25 of ischemia. Our results demonstrated 
that during myocardial ischemia phase Ia, the VW of a bidimensional tissue had an 
asymmetric unimodal distribution, peaking at minute eight after coronary occlusion (see 
Fig. 10D), which is in close agreement with experimental observations (Cascio, 2001). We 
also analyzed the SF for conduction in the regional ischemic tissue using an improved 
version (Romero et al., 2009) of the formulation proposed by Shaw and Rudy (Shaw & 
Rudy, 1997c). Our results also indicate that, as ischemia progresses, the SF decreased and 
the heterogeneity of refractoriness grows. Therefore, a lack of correlation between 
dispersion of refractoriness and vulnerability to reentry was observed, although our results 
corroborate the fact that dispersion of refractoriness is essential to reentry generation. 
Moreover, we found that, in approximately half of the simulated reentries, the area where 
unidirectional block (UDB) took place was completely recovered from refractoriness when 
the AP propagation failed. However, the line where the SF dropped below unity matched 
very much the area where propagation failed in every simulation. Thus, it seemed that the 
reduction of the source-sink ratio is the ultimate cause of the UDB leading to reentry. 
Reentrant patterns of activation during regional acute ischemia have been mechanistically 
analyzed also under the effects of pharmacologic agents, such as pinacidil or lidocaine. 
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Fig. 10. Panel A: Evolution of electrophysiological parameters after the onset of myocardial 
ischemia in the CZ: [K*]o, farr, and relative activation of INa (fna) and Icar (fca). Panel B: 
Regional ischemic tissue layout comprising a NZ, a BZ, and a CZ. Spatial gradients of 
electrophysiological parameters are described above the tissue. 


4.2 Drugs effects on reentry during myocardial ischemia 

To analyze the arrhythmogenic effects of pinacidil, different sets of simulations were carried 
out using different concentrations of the drug. In each set the regionally ischemic 2D tissue 
described above was stimulated using the aforementioned Sı-S2 protocol and the VW was 
evaluated. Table 1 shows that increasing the concentration of pinacidil had an interesting non- 
monotonic effect on the duration of the VW. Indeed, low concentrations of the drug increased 
the VW, until a maximum of 54 ms was reached for [P]=3 um/L. If the concentration was 
further increased, the VW decreased, until the VW vanished for [P]=10 um/L. 

The simulations allowed to unravel the mechanisms underlying pinacidil-induced changes in 
vulnerability, through the analysis of the patterns of excitation in the tissue and the time 
evolution of ionic currents and membrane potential. Our results showed the existence of two 
opposite effects of pinacidil in terms of AP propagation, which modulate the biphasic effects 
on the VW. On the one hand, pinacidil enhances Ixiarp) (Ferrero, Jr. et al, 1996), which 
counteracts inward currents responsible for depolarization, mainly Ina (which is depressed as 
a consequence of ischemia) and Ica), so that AP reaches lower values of Vm during the 
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depolarization phase. In this way, the potential gradient in the direction of propagation is also 
reduced, and less electrotonic current flows to downstream cells (Shaw & Rudy, 1997a), 
reducing axial current. In these conditions, the SF is reduced (Shaw & Rudy, 1997c) and 
propagation failure may develop. On the other hand, the Ikarp) enhancement provoked by 
pinacidil reduces APD, which results in an earlier recovery of excitability (Shaw & Rudy, 
1997b), favoring AP propagation. Both effects of pinacidil are related as follows: the earlier the 
cell recovers its excitability, the less axial current is needed to elicit an AP. 


Pinacidil (mmol/L) VW (ms) 
AT | sper ap) 
0,001 49 [166-215] 


0,003 54 [166-220] 
0,005 46 [179-225] 
0,007 14 [207-221] 


Table 1. Width of the vulnerable window (VW) for different concentrations of pinacidil. 
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Fig. 11. Width of the VW in ms for different degrees of ischemia (5, 6, and 7 minutes after 
occlusion) and concentrations of lidocaine (20, 50, and 100 umol/ L). 


We also simulated the effects of lidocaine on a regionally ischemic bidimensional tissue 
prone to reentrant circuits. Simulations were also carried out using an S1-S2 protocol applied 
on one edge of the tissue and the VW for reentry was calculated for various degrees of 
severity of ischemia and various concentrations of the drug (20, 50, and 100 umol/L). The 
VWs for reentry are shown in Fig. 11, and follow an unimodal behavior with ischemia time 
course, peaking (35 ms) for 6 minutes after the onset of ischemia. When the concentration of 
lidocaine was increased, the VW became wider, indicating a higher vulnerability to reentry 
in the presence of the drug. This effect can be due in part to the heterogeneous action of the 
drug in the diverse zones of the tissue, contributing to the dispersion of CV and setting the 
stage for reentry. 
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4.3 Effects of dofetilide on 1D heterogeneous tissue 

In the case of dofetilide, a class II antiarrhythmic drug recently approved by FDA (Food and 
Drug Administration) for the treatment of persistent atrial fibrillation and flutter, different 
studies have questioned the antiarrhythmic action of dofetilide in preventing and terminating 
ventricular tachycardias. Indeed, dofetilide promotes the prolongation of QT interval, which 
has been related to the trigger of a polymorphic ventricular tachycardia called torsade de 
pointes. The study of the proarrhythmic effects of dofetilide has also been undertaken with 
computer simulations at the tissue level, where pseudo ECGs have been evaluated, giving 
understanding about the mechanisms by which QT interval is altered by this drug. 
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Fig. 12. APDoo distribution in a heterogeneous linear strand, APs of three cells belonging to 
the different regions of the tissue and pseudoEGG, in two different conditions: before the 
application of dofetilide (control) and when the steady-state in drug-binding was reached 
after the application of dofetilide 100 nmol/L (BCL= 1000 ms), for normal coupled fiber 
(panel A) and under poor coupling conditions (panel B). 
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To study the effect of dofetilide on transmural dispersion of repolarization (TDR), a linear 
strand including the three kinds of cells was used. Fig. 12A illustrates the steady-state APDoo 
distribution along a well-coupled heterogeneous fiber and the steady-state APs time-course of 
the three types of cells belonging to the different regions of the tissue (right side in Fig. 12) 
under two different conditions: before (control) and after the application of dofetilide 100 
nmol/L. The pseudo ECGs are also represented before and after the application of the drug. 
Our results show that in a well-coupled strand under control conditions, at a BCL of 1000 
ms, M-cells present a maximum APDo9 of 191 ms, which is 15 ms and 35 ms longer than the 
shortest APD of endocardial (176 ms) and epicardial (156 ms) cells, respectively. The 
application of dofetilide 100 nmol/L increased APD dispersion (difference between the 
maximum and the minimum APDop in the fiber) from 35 ms to 58 ms, being the maximum 
APDo recorded in the M-zone and the minimum in the epicardial zone. The electrograms 
show an increment in the QT interval induced by the application of dofetilide from 220.8 ms 
(control) to 274.6 ms for 100 nmol/L dofetilide. 

The effect of decreasing the intercellular coupling (intercellular coupling resistance Rj 
multiplied by three) at a BCL of 1000 ms is shown in Fig. 12B. Poor coupling increased the 
APDo dispersion from 35 ms (normal coupling) to 58 ms. Even more, under poor coupling 
conditions dofetilide induced a steeper distribution of APD along the fiber and the APDoo 
dispersion increased to 102 ms, for dofetilide concentration of 100 nmol/L. The QT intervals 
were higher than the ones observed in a well coupled fiber. The application of 100 nmol/L 
dofetilide induced a sharper increase of QT from a value of 251.8 ms in control, to 318 ms. 


5. Modeling whole heart electrical activity 


Whole heart models are needed to study arrhythmias that critically depend on the spatial 
organization of the heart. Some spatial physiological features may be localization 
dependent, such as regional ischemia or localized infarct scars. Furthermore, certain 
arrhythmias are uniquely encountered in the whole heart, such as large reentrant circuits in 
acute ischemia or initiation of atrial fibrillation by pulmonary venous foci. Other important 
considerations of the whole heart models are the cardiac fibers arrangement and transmural 
inhomogeneity in repolarization. Cardiac fibers are arranged as counter-wound helices 
encircling the ventricular cavities, and the local orientation of these fibers depends on 
transmural localization (Helm et al., 2005). On the other hand, the heterogeneous expression 
of Ik, It, and Incx among others, leads to transmural inhomogeneity in repolarization 
(Antzelevitch & Fish, 2001; Antzelevitch et al., 1991). The vertiginous progress in medical 
images allows building realistic models of the heart including its highly complex anatomical 
structure of ventricles and/or atria (see Clayton et al. 2011 for review). 

The mathematical problem associated to the resolution of the differential equations 
modeling the electrical propagation in a 3D heart model does not have an analytical 
solution. Indeed, electric conduction is described by a set of partial differential equations 
(PDEs) and ionic currents through the cell mambrane are described with a nonlinear stiff 
system of ordinary differential equations (ODEs). The numerical solution of the equations is 
very computationally demanding (Heidenreich et al., 2010b) and require the discretization 
in space and time of PDEs, as well as the integration of nonlinear systems of ODEs. The 
operator splitting technique for the time discretization (Qu & Garfinkel, 1999) is used to 
solve the equations system. The numerical methods used, such as finite differences method 
or finite elements method, are computationally demanding and require the use of high 
performance computing techniques. 
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Whole heart models also permit to relate the arrhythmic behavior with its manifestation in 
electrograms on the cardiac surface and in electrograms on the body surface, through a torso 
model (Weiss et al., 2009). 


5.1 3D simulations of myocardial ischemia 

The use of these ventricular 3D models has shed light into the mechanisms by which 
reentrant arrhythmias are initiated during regional myocardial ischemia. As an example, a 
3D model of the human ventricles has been recently used to study the appearance of figure- 
of-eight reentry under regional acute ischemic conditions (Heidenreich et al., 2010a; 
Heidenreich et al., 2010b). 

Fig. 13 shows different views of the simulated human ventricles. The detailed realistic local 
fiber orientation obtained from DTI images (diffusion tensor imaging) was introduced in the 
model and is key in determining the electrical propagation. An ischemic region was 
introduced in the left ventricle, following the same normal-border-central zone scheme 
explained in a previous section of this chapter. The central and right panels in Fig. 13 depict 
the values and gradients of the main ischemic parameters (namely [K+]o [ATP], Ina 
maximum conductance and Ica) maximum conductance). 
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Fig. 13. Anatomically realistic human ventricles subject to regional acutely ischemic 
conditions (Heidenreich et al., 2010b). Top A panel: render image of the simulated 
ventricles. Bottom A panel: fiber orientation as obtained from DTI. Panel B: distribution of 
[K+]o [ATP]i, Ina and Icaq) in the CZ, the NZ and the different epicardial BZs. 
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Fig. 14. Electrical activity in the simulated 3D human ventricles shown in Fig. 13 after the 
delivery of a premature extrastimuls in the BZ after a train of regular conditioning stimuli. 
A figure-of-eight reentry is initiated, which is very similar to experimentalty observed 
ischemic reentry (Janse & Wit, 1989). 


As shown in Fig. 14, the model predicted the generation of figure-of-eight reentries, which 
cross the central ischemic zone formed in the epicardial surface due to the longer refractory 
period of the midmyocardial layers. Also, focal activity experimentally observed in the 
epicardium could be caused by re-entrant wavefronts propagating in the mid-myocardium 
that reemerge in the heart surface. 


6. Conclusion 


Multiscale modeling of the electrical cardiac activity represents nowadays a very valuable 
tool in cardiac disease research which complements experimental and clinical studies. 
Although significant progresses are made in this field, further improvements are required, 
such as the modeling of electromechanical coupling of the heart activity, which currently 
represents an important challenge for modelers. The vertiginous technical improvements in 
the obtention of medical images and genetic, molecular and ionic measurements, represent a 
great help to the modeling field, which is in the process of creating a virtual human. 
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